# Figure 1
rm(list=ls())
gc()

quartz(width = 10, height = 5, type = "pdf", "Figure0", file = "Figure0.pdf")

par (fig=c(0,1,0,1), # Figure region in the device display region (x1,x2,y1,y2)
	omi=c(.8,0,0,.5), # global margins in inches (bottom, left, top, right)
	mai=c(0,.8,1,.3)) # subplot margins in inches (bottom, left, top, right)
	layout(matrix(c(1,2), 1, 2, byrow = TRUE), widths=c(3,3), heights=c(3,3))

## Panel 1: Historical and forecast temperatures
hadcrut <- read.table("HadCRUT.4.5.0.0.annual_ns_avg.txt")
hadcrut <- hadcrut[,1:2]

cmip <- read.csv("das_CMIP5_rcp8_5.txt", stringsAsFactors=F)

plot(hadcrut[,1], hadcrut[,2]+1, type="l", ylim=c(0,10),xlim=c(1850,2100), xaxs="i",yaxs="i",bty="L",yaxt="n",xlab="",ylab="")
axis(side=2,at=seq(0,10,2),las=2)
lines(cmip[2:nrow(cmip),1], cmip[2:nrow(cmip),42]+2)
names(cmip)[42]
cmip[1,42]

lines(cmip[2:nrow(cmip),1], cmip[2:nrow(cmip),2]-hadcrut[1,2]+3)
names(cmip)[2]
cmip[1,2]

lines(cmip[2:nrow(cmip),1], cmip[2:nrow(cmip),22]-hadcrut[1,2]+4)
names(cmip)[22]
cmip[1,22]

mtext(side=1, "Year", outer=T, line=2, cex=.9, at=0.275)
mtext(side=3, expression(Delta * "T + offset (˚C)" ), outer=T, line=-4, cex=.9, at=0.07)

## Temperature distributions
mu = 8
sd = .4

t <- seq(0,100,.01)
t1 <- cbind(t,dnorm(t, mean = 2, sd = .1))
t2 <- cbind(t,dnorm(t, mean = 2, sd = sd))
t3 <- cbind(t,dnorm(t, mean = mu, sd = .1))
t4 <- cbind(t,dnorm(t, mean = mu, sd = sd))


# Panel 2: Reflecting temperature distribution through damage function
# Damage distributions
d4 <- 10*(1-(1/(1+(t/20.46)^2 + (t/6.081)^6)))
d4 <- cbind(d4, dnorm(t, mean = mu, sd = sd))

plot(t1[,1],t1[,2],type="n",xlim=c(0,10),ylim=c(0,10),xaxs="i",yaxs="i",bty="L",yaxt="n",xaxt="n",xlab="",ylab="")
axis(side=1,at=seq(0,10,2),label=seq(0,1,.2))
axis(side=2,at=seq(0,10,2),las=2)
mtext(side=3, expression(Delta * "T (˚C)"), outer=T, line=-4, cex=.9, at=0.55)

mtext(side=1, expression("Damages as share of output, D(" * Delta * "T)"), outer=T, line=2, cex=.9, at=0.8)

tmin <- mu - 3*sd
tmax <- mu + 3*sd
dmax <- 10*(1-(1/(1+(tmax/20.46)^2 + (tmax/6.081)^6)))
dmin <- 10*(1-(1/(1+(tmin/20.46)^2 + (tmin/6.081)^6)))

polygon(x=c(0,dmax,dmax,0),y=c(tmin,tmin,tmax,tmax),col="grey85",lty=0)
polygon(x=c(dmin,dmin,dmax,dmax),y=c(0,tmax,tmax,0),col="grey85",lty=0)

tmin <- mu - 2*sd
tmax <- mu + 2*sd
dmax <- 10*(1-(1/(1+(tmax/20.46)^2 + (tmax/6.081)^6)))
dmin <- 10*(1-(1/(1+(tmin/20.46)^2 + (tmin/6.081)^6)))

polygon(x=c(0,dmax,dmax,0),y=c(tmin,tmin,tmax,tmax),col="grey",lty=0)
polygon(x=c(dmin,dmin,dmax,dmax),y=c(0,tmax,tmax,0),col="grey",lty=0)

tmin <- mu - sd
tmax <- mu + sd
dmax <- 10*(1-(1/(1+(tmax/20.46)^2 + (tmax/6.081)^6)))
dmin <- 10*(1-(1/(1+(tmin/20.46)^2 + (tmin/6.081)^6)))

polygon(x=c(0,dmax,dmax,0),y=c(tmin,tmin,tmax,tmax),col="grey60",lty=0)
polygon(x=c(dmin,dmin,dmax,dmax),y=c(0,tmax,tmax,0),col="grey60",lty=0)

## Temperature distribution
lines(t4[,2],t4[,1])

## Damage function
lines(10*(1-(1/(1+(t/20.46)^2 + (t/6.081)^6))),t) # Weitzman's damage function

## Damage distribution
lines(d4[,1],d4[,2])

dev.off()